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ABSTRACT 

The discovery of over 50 planets around evolved stars and more than 35 debris discs orbiting 
white dwarfs highlight the increasing need to understand small body evolution around both 
early and asymptotic giant branch (GB) stars. Pebbles and asteroids are susceptible to strong 
accelerations from the intense luminosity and winds of GB stars. Here, we establish equations 
that can model time-varying GB stellar radiation, wind drag and mass loss. We derive the 
complete three-dimensional equations of motion in orbital elements due to (1) the Epstein and 
Stokes regimes of stellar wind drag, (2) Poynting-Robertson drag, and (3) the Yarkovsky drift 
with seasonal and diurnal components. We prove through averaging that the potential secular 
eccentricity and inclination excitation due to Yarkovsky drift can exceed that from Poynting- 
Robertson drag and radiation pressure by at least three orders of magnitude, possibly flinging 
asteroids which survive YORP spin-up into a widely dispersed cloud around the resulting 
white dwarf. The GB Yarkovsky effect alone may change an asteroid’s orbital eccentricity 
by ten per cent in just one Myr. Damping perturbations from stellar wind drag can be just 
as extreme, but are strongly dependent on the highly uncertain local gas density and mean 
free path length. We conclude that GB radiative and wind effects must be considered when 
modelling the post-main-sequence evolution of bodies smaller than about 1000 km. 

Key words: minor planets, asteroids: general - Kuiper belt: general stars: AGB and post- 
AGB - stars: evolution - stars: white dwarfs - planets and satellites: dynamical evolution 
and stability 


1 INTRODUCTION 

Robust observational evidence for exoplanetary systems around 
both main sequence (MS) and post-MS stars demands a self- 
consistent theory for the formation and late evolution of plane¬ 
tary systems. Although the predominant focus of exoplanetary 
science has been MS systems, the last decade has seen a surge 
in interest for giant branch (GB) stars and white dwarfs (WD) 
that host planetary systems. 


1.1 Importance of post-MS planetary systems 

The discovery of over 50 exoplanets and nearly 100 substel- 
lar companions to stars which have left the main sequence (see 
IWang et al.ll2014l . and references therein) help motivate studies 
which attempt to link the past and future evolution of plan¬ 
etary systems. Detections of these companions have become 
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more robust with the combined efforts of Dopp ler radial veloc¬ 
ity observations and transit-based photom etry dLillo-Box et al.l 
120141 ; ICiceri et al JI2015I : [Ortiz_et i ^ih 2015k_Some GB stars con¬ 
tain multiple planets (e.g. Niedzielski et al.ll2015l l. and/or en¬ 
hanced abundances of lithi um, which could represent a key 
tracer of recent engulfmen t dAdamow et al.ll2012l : iNowak et al.l 
120131 : lAdamow et al.l[20l4). Observations of dust y debris discs 
around subgiant stars ( Bonsor et al.l 120131 . 120141) suggest that 
collections of material, such as the asteroid belt, can survive the 
entire MS lifetime of a star. 

These observations are bolstered by strong supporting 
evidence for asteroidal material in WD systems. Although 
only a few dusty debris discs are known around subgiant 
stars, over 35 of such discs have been observed o r biting WDs 
|Zu^ffirman_&^ecklin [1987Becklimetah|_2005kKihcet_aH 


20051: iReach et al.ll2005l: iFarihi et al.ll2009l: iBergfors et al.ll2014l: 


gaseous components 


Rocchetto et al] 201511. Sometimes accompanying the dust are 
lIGansicke et ai1l2006l. 20071 20081: Gansickel 


1201 ll ; iDufour et al.l [20121 : IFarihi et al.ll2012l ; |Melis et al.| |2012|) 

with both components overlapping radially ( Brinkworth et ahl 
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cally does not exceed one Solar radius, which is close to the 
tidal disruption radius for an asteroid. The environment within 
this compact region of space is dynamic: the gaseous compo¬ 
nents of the discs provide kinematic information, and demon- 
strate sometimes extr e me variations with e very new observatio n 
llGansicke et al.l 120081 : IWilson et al.l 120141 ; iManser et al.l l2015l l . 
Sharp decreases in luminosity in dusty disc systems also indi¬ 
cate qu ick drastic change s in the morphology of the circumstellar 
debris (IXu fe Jurall20l4) . 

WDs accrete the disc material, which appear as atmo¬ 
spheric metal abundances, or “pollution”. Any metals heavier 
than helium sink out of sight on time scales of days to Myrs 
(depending on the depth of the convection zone), which are in 
all cases much shorter than the WD cooling age, i .e the tim e 
since the star became a WD (see Fig. 1 of lWvatt et alj|2014h . 
Therefore, The detection of photospheric metals implies recent, 
or ongoing accretion from a circumstellar debris reservoir. Fur¬ 
ther, between 25 to 50 per cent of all WDs con t ain signatures 
of me tal pollution (IZuckerman et aT]|2003l . l20ld : iKoester et al.l 
l20l4l . and at widely varying (over 5 Gyr) cooling ages. 


1.2 The link between GB and WD planetary systems 


These ubiquitously active environments close to WDs must be 
linked with planetary system evolution during the GB phases 
of evolution. The interstellar medium is too diffuse to ac¬ 
count for amount of metals seen, and is predominately com¬ 
posed of hydrogen. Atmospheric abundances of WDs which 
are hydrogen-poor and metal-rich (known as DBZ WDs) can¬ 
not then arise from the in t erstel l ar medium jA anncstad et_ah| 
19931: [Friedrich et al.l 120041: I Jura 120061 : iKilic fe Rcdficld 20071 : 


Farihi et al.1 20ld : Barstow et akf 20141 ') . 


The currently-favoured explanation for the origin of atmo¬ 
spheric pollution and debris discs are asteroids which are dy¬ 
namically perturbed close to the WD and then tidally disrupted 
there. By asteroids we refer to any objects whose radius is ap¬ 
proximately between O.lkm-lOOOkm lying within thousands of 
au of their parent stars; we denote pebbles and planets as objects 
with radii between 1mm- lm and larger than 1000km, respec¬ 
tively. Evidence suggests that other potential explanations are 
unlikely. The accretion and disruption of planets around WDs 
likely does occur, but too infrequently to explain th e large frac¬ 
tion o f all WDs currently-observed to be accreting (IVeras et alj 
l2013al : iMustill et al.l 12014 IVeras fe Gansickj l2015h . Exo-Oort 
cloud comets (with semimajor axes of tens of thousands or hun¬ 
dreds of thous ands of au) are unlike ly debris progenitors on 
both chemical iTZuckerman et al.ll2007l l and dynamical grounds 
llVeras et al.ll2014al ), altho ugh in isolated ca ses they may pro¬ 
duce detectable accretion fStone et al.ll2015r ). 

Asteroids provide a readily-available reservoir of material 
with a range of metals diverse enough to help ex plain the ele¬ 
mental medley observed in a tmospheric poll ution (iDufour et al.l 
120121 : iGansicke et al.l 120121 : IXu et al] l20l4) . Further, the wa¬ 
ter rete ntion level in asteroids during the GB phases is un¬ 
certain (lJura fe Xul l2010l . 120121 ') , but might be high enough to 
explain WDs wit h water-rich atmospheres (IFarihi et al.ll2013l : 
iRaddi et al.ll2015l '). 

Dynamically, asteroids can be scattered onto or close to the 


WD through gravitational interactions with surviving planets. 
Subsequently, th e asteroids wh i ch av o id direct co l lisions might 


tidally disrupt (Graham 

et al. 199C 

: Jura 

120031: Debes et al.1 

20lJ; Bear & Sokel20lJ 

Veras et al. 

2014b 

), creating discs and 

accreting onto the WD ( 

Rahkov 2011a.hi: 

Rahkov & Garmillal 


120121 ; IWvatt et al.l I20l4) . Some investigations which modelled 
asteroid perturbations during and beyond GB mass loss have 
explored different regions of the available phase space, and in- 
clud ed one planet in th eir simulations. These studies encom¬ 
pass iBonsor et al.l d201 ll j , who considered asteroids in an exo - 
Kui per belt (at 30 au from the star), and iDebes et alj (120121 s ) 
and iFrewen fe Hansenl l|20l4) . who modelled a belt more akin 
to the asteroid belt (at about a few au from the star). 


1.3 Previous most relevant work 

None of those three studies included effects from radiation or 
stellar wind drag in their numerical simulations during GB evo¬ 
lution. However, IVeras et al.l (12014 jj demonstrated that nearly 
all asteroids with radii between 100m and 10km within about 7 
au of their parent star on the MS will later be destroyed because 
of rotational fission - due to the GB star’s luminosity - from the 
YORP effect. Consequently, an exo-belt similar to the asteroid 
belt simply will not survive the GB stellar phases intact, and 
this phenomenon should be taken into account when modelling 
asteroid evolution during the WD phases. 

Yet, YORP-induced fission is just one c onsequence of the 
violen t dynamical environment of GB stars. IBonsor fe Wvatt] 
(120101 s ) described several others in the context of their influ¬ 
ence on debris discs orbiting stars which turn off of the main 
sequence. The authors considered the effects from Poynting- 
Robertson drag, and forces they denoted as radiation pressure, 
stellar wind pressure, and stellar wind drag. Here, we isolate 
these forces and derive how they change the orbits of small bod¬ 
ies in a general fashion, without making any assumptions about 
size distributions other than that their diameters are large com¬ 
pared to the wavelength of the incident radiation. We determine, 
for example , the eccentricity a nd inc lination evolution of orbit¬ 
ing bodies; IBonsor fe Wvattl (l2010l l assumed that the bodies 
were on circular co planar orbits. 

In addition to IBonsor fe Wvattl (|2010T ). Ibong et all (120101 s ) 
also considered physical effects arising from GB stars. 
iDong et ~aL 1 d2010l) focused on stellar wind drag and entrainment 
of small particles along the GB phases, and explicitly included 
drag in their equations of motion (Stokes regime only), but did 
not present those equations in orbital elements. They neglected 
Poynting-Robertson drag, YORP and Yarkovsk^j effects during 
AGB evolution, but estimated that those effects would lead to 
significant orbital evolution - a notion that we quantify in detail 
here. 

Before proceeding with the body of the paper, we present a 
summary chart (Fig. [TJ of the physical effects which should be 
considered when modeling pebbles, asteroids or planets along 
different phases of stellar evolution. Boxes without checkmarks 
indicate forces that are negligible along that particular phase. 

This chart simultaneously showcases our conclusions and 


1 IDong et al.l J2010I ') do not make a distinction between the Yarkovsky 
and YORP effects. 
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Figure 1. Effects which are necessary to consider for modeling. The 
designations MS, GB and WD refer to the main sequence, giant 
branch and white dwarf phases of stellar evolution. The radiative 
forces include those from “PR-drag and RP” (Poynting-Robertson 
drag and radiation pressure; see Section 2.5), the Yarkovsky effect 
(see Section 2.6), and the YORP effect (see Section 2.2). Properties 
of stellar winds are introduced in Section 3.3. 


reaffirms some re sults by both iBonsor fc Wvattl d20ld ) and 
iDong et al.l d20ld ). For example, both studies find that the dis¬ 
tribution of pebbles in WD systems will be crucially determined 
by wind drag along the GB phase. Multiple checkmarks in a 
single row of Fig. |T] indicate either that multiple forces must 
be considered together for that stellar phase, or that one may 
dominate the motion, depending on the physical and orbital 
variables chosen. Additional forces not in cluded in the chart, 
such as collisio ns dBonsor fe Wvattl l20ld) or those from addi¬ 
tional planets (|Dong et al.ll20inl ). can alter the dominant force 
in a given system during a given phase. 


1.4 Outline for this article 

Here we determine how GB stars change the orbits of bodies 
in an all-inclusive fashion under the guise of the perturbed two- 
body problem. We perform direct comparisons of the effects due 
to gravity, radiation and drag. We then derive equations of mo¬ 
tion in orbital e lements throug h the eq uation-generation proce¬ 
dure outlined in lVeras fe Evana (l2013al ). All of the variables and 
parameters used are delineated in Tables 1-2 for easy reference. 
Henceforth, we denote the bodies as targets to help avoid any 
bias towards pebbles, asteroids or planets. Let m represent the 
target’s fixed mass, v the velocity of the target with respect to 
the star, and r the distance between the centre of the target and 


the centre of the stafl The star has a mass of M which does 
change with time. The target then evolves according to 

dv _ G(M(t = 0) + m)r 
dt r 3 



where the three rightmost terms of equation © refer to the 
contributions from stellar mass loss, radiation, and wind drag. 
Cartesian formulae for wind drag and mass loss are well- 
established; radiative effects provide more subtleties. Complete 
sets of orbital element formulae have been previously described 
only for mass loss, to our knowledge. 

In Section 2, we derive the nontrivial expression for 
(du/df) ra in terms of Cartesian positions and velocities by de¬ 
veloping a framework in which to consistently treat Poynting- 
Robertson drag, radiation pressure, the Yarkovsky effect and 
YORP. Previous investigations almost exclusively focused on 
the Solar system. Our treatment here removes the biases of con¬ 
stant luminosity and asteroid belt location which are inherent 
in those studies. 

Section 3 contains our comparison of the magnitudes and 
relevant size regimes of all of the terms in equation ©. These 
comparisons of instantaneous accelerations help quantify the rel¬ 
ative importance of the different forces at a given time, but not 
over long times (secular timescales). 

In Section 4 we derive the resulting (unaveraged) equations 
of motion in orbital elements for all of these effects, when viable. 
We then derive the averaged equations of motion, which reveal 
secular aspects of the motion. These equations provide a way 
for us to either compute specific trajectories or place bounds 
on the motion given our framework in Section 2. We summarize 
our findings in Section 5. 


2 A FRAMEWORK TO MODEL RADIATIVE 
EFFECTS 

Over the last century, the movement of a Solar system ob¬ 
ject due to radiation pressure has popularly bee n described 
throu gh Poynting-Robertson drag rtPovntinef 1904; Robertsonl 
Il937l ). the Yarkovsky effect llRadzievskiil 1 1954 : IPetersonl Il97^ ) 
and the YORP (Yarkovsky-O’Kee f e-Radzievskii-P addack) effect 
dRadzievskiil 1 1 9541 : lO’Keefd Il976l : IPaddackl 1 19691 ) @. Poynting- 
Robertson drag and the Yarko vsky effect, whi ch are rarely men¬ 
tioned in the same context (see |Petersonlll976l for an exception), 


2 This distance is between the centres of the bodies and not their 
surfaces because we will assume Gauss’ law to calculate the incident 
radiation on the asteroid. Therefore, as long as the target remains 
outside of the star’s surface and the stellar emission does not signifi¬ 
cantly vary over the star’s surface, our equations will hold. Further, 
the momentum transport equations contain the equivalent cross sec¬ 
tion of the asteroid, which is assumed to be located where the asteroid 
has its largest diameter. For a spherical asteroid the cross sectional 
plane contains the asteroid’s centre. 

3 Yarkovsky’s original pap er, dating from around 1900, has appar¬ 
ently been lost dQpiklll95lh ■ 
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Table 1. Roman variables and parameters used in this paper. 


Variable 

Explanation 

Equation # 

a 

Target’s semimajor axis 


A 

Target’s cross-sectional area (assumed to be = i tR 2 here) 


B 

Auxiliary piecewise variable 

ESI 

c 

Speed of light 


Ci - Cn 

Auxiliary variables 

HOES] 

c 

Target’s specific heat capacity 


d 

Target’s diameter 


D 

Auxiliary vectors 

1861921 

e 

Target’s eccentricity 


f 

Target’s true anomaly 


V 

Target’s penetration depth scaling 

El 

G 

Gravitational constant 


h 

Target’s specific orbital angular momentum (= r x v) 


i 

Target’s inclination with respect to the star’s equator 


n 

3x3 unit matrix 


k 

Constant between 0 and 1/4 

EE 

K 

Target’s thermal conductivity 


K 

Constant equal to 0.165 


L 

Star’s luminosity 


M 

Star’s mass 


m 

Target’s mass 


mu 

Mass of Hydrogen atom (fx 1.66 X 10 -27 kg) 


n 

Target’s mean motion 


p 

Momentum 


p 

Power 


Qahs 

Target’s absorption efficiency (dimensionless) 


Q ref 

Target’s reflecting efficiency, or albedo (dimensionless) 


Q yar 

Difference of Q a b s an d Qref (dimensionless) 


Qpr 

Target’s radiative efficiency due to Poynting-Robertson drag (dimensionless) 

El 

Q 

3x3 diagonal radiation matrix 

EE 

r 

Distance between the centre of the target and the centre of the star 


R 

Target’s radius 


R* 

Star’s radius 


n 

Auxiliary variable 

El 

R 

General 3x3 rotational orbital matrix 

1401481 

Ri (s) 

3x3 Rotational matrix #1 for target spin 

nsi 

R 2 (s) 

3x3 Rotational matrix #2 for target spin 

EE 

Ri (h) 

3x3 Rotational matrix #1 for orbital angular momentum 

EE 

R 2 (h) 

3x3 Rotational matrix #2 for orbital angular momentum 

EE 

My (s* 5 </>) 

3x3 Rotational matrix for total diurnal Yarkovsky contribution 

EE 

r y (h,A 

3x3 Rotational matrix for total seasonal Yarkovsky contribution 

EE 

Re 

Reynolds number of stellar wind at the target’s location 

EE 

s 

Target’s spin axis (with angular speed = |sj) 

El 

^crit 

Angular speed beyond which bodies in the Solar system with 0.1 km < R < 10.0 break up 

El 

T 

Temperature 

EE 

V o 

A fiducial Solar force (= 10 17 kg m/s 2 ) 

El 

V 

Target’s velocity with respect to the star 


% 

Gas velocity (where the gas is the stellar wind) 


V s 

Local sound speed of gas (stellar wind) 


Vrot 

Star’s rotational velocity at its equator 


w 

Coefficient of the target’s physical asymmetry (dimensionless) 

El 

Y 

3x3 diagonal Yarkovsky matrix 

EE 

Z 

Stellar latitude 

© fgp RA 
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Table 2. Greek variables and parameters used in this paper. 


Variable 

Explanation 

Equation # 

p 

Ratio of target’s acceleration from radiation to that from gravity 

ED 

7 

Ratio of target’s acceleration from radiation to that from mass loss 

m 

<5 

Ratio of target’s acceleration from radiation to that from wind drag 

1571581 

e 

Target’s emissivity 


c 

Mean free path length of gas (stellar wind) 


e 

Collision cross-section of gas molecules 


t 

Relativistically-corrected direction of incoming radiation 

0 

K 

Auxiliary variable 

[All 

A 

The peak wavelength of absorbed radiation 


P 

Mean molecular weight 


5 

Seasonal thermal lag angle in the orbital plane 


n 

Target’s orbital period (= 2tt/ n) 


p 

Target’s mass density 


Pg 

Mass density of gas (stellar wind) 


a 

The Stefan-Boltzmann constant 


s 

Target’s rotational spin period (= 2n/\s\) 


Ti -Tg 

Auxiliary variables 

lA6l 

<t> 

Dirurnal thermal lag angle along the target’s equator 


4- 

Mass loss index [= ( dM/dt)/(nM ) ] 

EH 

UJ 

Target’s argument of pericentre 


n 

Target’s longitude of ascending node 



both refer to a type of orbital recoil effect. YORP instead de¬ 
scribes the change in the spin of a body due to absorbed radia¬ 
tion. These notions are quantified more precisely in the following 
subsections. 

Although all three effects were originally defined in a Solar 
system context, their usage and names have been extended to 
extrasolar systems, and pa rt icularly post-main- sequence exosys¬ 
tems fe.g. I Rafikovl 1201 1 al ibi: IVeras et aI1l2014ch . Their applica¬ 
bility to extrasolar systems is at a fundamentally different level 


partly because of the difference in the 
Numerous studies (e.e. Krvszczyriska 

precision of available data. 
20131: Rozitis et all 20131: 

Jacobson et al.l 

2014l;lLowrv et al.l 201 

^Lupishko & Tielieusoval 

20 lT Polishook 

2014: Bottke et alj 20 

15f) performed simulations 


to match outcomes to known solar system objects at a standard 
which will be unobtainable in exosystems in the foreseeable fu¬ 
ture. Therefore, this paper emphasizes placing bounds on the 
motion and learning as much as possible from the equations of 
motion, rather than concerning itself with detailed simulation 
suites. But first, we must create a self-consistent model that is 
free from the Solar system biases of a constant-luminosity star 
and asteroids concentrated in location-specific families. 


2.1 Goal 

Our goal in Section 2 is to find (du/d£) ra , which can be divided 
into the following three terms 



The first term on the RHS of equation @ refers to the extra ac¬ 
celeration on the target due to the absorption of radiation. The 


second and third terms together give the acceleration due to the 
re-emission of radiation. The second term refers to re-emission 
from immediate reflection, and the third term to re-emission 
from delayed reflection. This delayed reflection is the Yarkovsky 
effect, and arises from internal thermal redistribution of the ab¬ 
sorbed radiation. Note that YORP is absent from equation m 
because it provides no direct orbital acceleration to the target. 


2.2 Acceleration due to YORP 


Nevertheless, crucially, YORP spins up or down the target. As¬ 
sume that the target rotates at an angular speed s around its 
spin axis, defined as s = ts x ,s y ,s z | T . Then Y ORP causes the 
angular speed to evolve according to dScheeresll2007i l 


= w ( i A (rr m\ 

dt 2npR 2 V ® L& ) 


(3) 


where we gave the star’s luminosity L a time dependence for 
emphasis, 0 ^ W < 1 is a constant determined by the physical 
properties (roughly the asymmetry) of the target, R is the tar¬ 
get’s radius, p is the target’s density, a is the semimajor axis of 
the orbit, and e is the eccen tricity of the orbit. The term in the 
last parenthesis comes from IVeras et all (l2014d 'l and represents 
a fiducial Solar force ( Uq = 10 17 kg m/s 2 ) for application to 
general extrasolar systems. 

Equation ([3j is important because it helps determine if a 
target will undergo fission or survive. In the latter case, the 
spin evolution could affect the target’s acceleration through the 
Yarkovsky effect, as described in the below sections. Solar sys¬ 
tem observations show a sharp cutoff for the maximum spin of 
asteroids with radii approximately between 0.1 km and 10 km 
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llHarrij|l994l : ljacobson et al.| [2014l. This critical spin, s cr it, cor¬ 
responds to a critical breakup period E cr it = 27r/s cr it ss 2.33 
hrs. If the target is assumed to be a strengthless rubble pilfl 
then we can express s cr it in terms of target density p as 


2tt 


^crit — 


V3n/Gp 


= 7.48 x 10" 


2 g/cn 


2 rad 
s 


(4) 


The interplay between the YORP and Yarkovsky effects 
is complex. For example, YORP may spin up a target quickly 
enough to activate the Yarkovsky effect. However, if the target 
is spinning too quickly, it will lose its thermal gradient (or break 
up). We do not consider the time evolution of the spin in this 
work. 


2.3 Acceleration due to absorbed radiation 

The expression for t he accelerat ion due to absorbed radiation is 
from equation 2 of I Burns et al . I (EH) 

/<Zv\ _ AL(t)Q ab3 _ . , 

\dt ) abs Anmcr 2 

where A is the geometric cross-sectional area of the illuminated 
surface, Q a b s represents the target’s absorption efficiency (or, 
equivalently, fractional amount of energy absorbed that con¬ 
tributes to the momentum transfer), c is the speed of light, and 
T is the relativistically-corrected direction of the incoming radi¬ 
ation to the target such that 



Our equation © differs from equation 2 from iBurns et al.l 
(|2014|) because we consider a general extrasolar system with a 
time-dependent lu minosity instead of the Solar constant. Also, 
IBurns et al.l |2QlII) explain that the primed quantities used in 
their equation are equivalent to the unprimed values in equa¬ 
tion © to low order in v /c ; this difference in quantities has 
previously been identified dKlacka et al . 1 EH). 


2.4 Acceleration due to immediately-reflected 
radiation 

The acceleration on the target due to immediately-reflected ra¬ 
diation is expressed in a similar form 

/ di A _ AL(t)Q ie { , 7 > 

\dt) ref 4nmcr 2 

where Q re f is the target’s reflecting efficiency, or, the fractional 
amount of energy immediately reemitted after absorption, or 
simply the albedo. In order to avoid having to account for com¬ 
plex scattering properties of tiny particles, we assume that our 
targets are at least 100 times the wavelength of radiation (cor¬ 
responding to the approximate limit between geometric optics 
and Mie scattering regime). 


2.5 Acceleration due to Poynting-Robertson drag 
and “radiation pressure” 

IBurns et al. idH) observed that the terms Poynting-Robertson 
drag and radiation pressure are ambiguous. Despite how the 
pressure created by radiation is responsible for all of the terms 
on the RHS of equation ©, plus YORP, the term radiation 
pressure has typically come to mean just the radial component 
of the sum of ()j|) abs and (^) re f • Similarly, the term Poynting- 
Robertson drag has typically come to mean just the tangential 
component of the sum of ( 7 ^) and (4 rr) .■ Thes e compo¬ 
nents are is olated and discu ssed in lBurns et all (|l979l l , iMienardl 
dl984l l and iHamiltonl dl993h . 

Also identified in previous studies is a Q coefficient that is 
associated with Poynting-Robertson drag, and is already aver- 
aged over the ste llar spectrum (see e.g. equations 4 and 19 of 
IBurns et al. IHzi). This coefficient is related to our formalism 
through: 

QpR — Qabs T Q ref- (8) 

One must keep in mind that the values of Q are averaged over 
the stellar spectrum because the stellar spectrum changes with 
stellar evolutionary phase. Such changes will not cause our con¬ 
clusions to deviate significantly. 


2.6 Acceleration due to the Yarkovsky effect 

Now we come to describing the last piece of equation ©, the 
perturbative acceleration due to the Yarkovsky effect. Recall 
that the Yarkovsky effect is due to reemission from delayed re¬ 
flection. The delay comes from the redistribution of thermal 
energy within the target. Therefore, if the target is too small, 
or spins very rapidly, no significant temperature differences will 
occur, and the Yarkovsky effect “turns off”. 


2.6.1 Penetration depth 


The minimum size of the target at which the effect becomes neg¬ 
ligible is approximately the penetration depth scale of the radi¬ 
ation. This depth scale can be derived from a he at conductio n 
model. Here we adopt the simplified ID model of Bro 1 (1200614 . 
This model provides the time lag between insolation and re¬ 
emission in realistic materials. The model also assumes that the 
target rotates at a constant angular speed, and therefore strictly 
equation © cannot be used in conjunction. Realistically, how¬ 
ever, YORP spin changes occur relatively slowly compared to 
thermal wave properties. Therefore, the YORP effect may be 
included in the model if s'is considered to be a function of time. 

This penetration depth scaling, V, is 



(9) 


where K is the thermal conductivity and C is the specific thermal 
capacity (or specific heat capacity). Hence, the Yarkovsky effect 
will “turn on” if the target diameter is larger thar0 V. As this 


4 Recently, I Sanchez fc Scheeresl (12014) derived spin barriers for bod¬ 
ies which harbour weak cohesive strength due to van der Waals forces 
between constituent grains. The critical breakup period is similar, but 
not monotonic as a function of decreasing target radius. 


5 In principle the thermal wave can heat the target by reaching down 
to its core. In this case, the sufficient condition for the Yarkovsky 
effect to “turn on” is if the target radius is larger than T>. 
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Figure 2. Minimum diameter of target at which the Yarkovsky effect “turns on”. This critical diameter is taken to equal the penetration depth 
scaling T>, which is a function of the thermal conductivity K and thermal capacity C. The left panel describes a high density p target with a fast 
spin period, equal to 1. 5E cr ;i. The right panel instead describes a low density slowly spinning target with E = 4.0E cr it. The plots bound realistic 
values of T> to between about 1 cm and 10 m. 



Figure 3. Relevant angles and vectors for the target. The spin axis 
is s, the specific angular momentum axis is h. 0 is the thermal lag 
angle along the target’s equator, and £ is the thermal lag angle in the 
orbital plane. The white spot corresponds to the direction of the star 
(subsolar point), and the red spot to the maximum of the thermal 
wave (the direction opposite to the Yarkovsky drift). 


scaling significantly depends on the target’s physical properties, 
we provide contour plots in Fig. 0 The contour plots cover the 
entire range of possible T> values for both a high density, quickly 
spinning (left panel) target and a low density, slowly spinning 
target (right panel). Note that 8 g/cm 3 is slightly higher than 
the density of iron and 0.2 g/cm 3 is slightly less dense than 
cork. Target thermal conductivities are confined to 10~ 2 — 10 3 
W/(mxK), which include all substances from air to silver. Tar¬ 
get specific heat capacities are confined to about 100 J/(kgxK) 
for heavy metals to about 2000 J/(kgxK) for water. The result, 
in all cas es, is 1 cm < T) < 10 m, conforming well to Figs. 25-26 
of I ro 1 (hood) . Therefore, we do not expect Yarkovsky to be 
active in sub-cm-sized targets, and we expect Yarkovsky to be 
always active in moderately fast spinning targets larger than 10 
m. 


2.6.2 Yarkovsky direction 


The geometry of the thermal redistribution determines the di¬ 
rection of the Yarkovsky acceleration. This redistribution is de¬ 
pendent on the orientations of both the spin axis s and the 
specific orbital angular momentum axis h = r x v. Let </> and £ 
represent lag angles of the thermal wave. The variable </ repre¬ 
sents the thermal lag angle along the target’s equator, i.e. how 
far the heat wave is dragged along from the subsolar point due 
to the target’s rotation. The variable £ represents the thermal 
lag angle in the orbital plane measured from the subsolar point. 
See Fig. [3] for a visual representation of these angles. 

Solar system studies traditionally denoted the contribution 
from the spin axis and 4> as “diurnal” and the contribution from 
the specific angular momentum axis and £ as “seasonal”. These 
characterizations can be traced historically to the thermal im¬ 
balance between the morning and afternoon hemispheres (for 
diurnal) or the winter and summer hemispheres (for seasonal) 
of a Solar system asteroid. Mathematically, the contributions 
are expressed via rotation matricies through an angle (here (j) 
and 0 about a given axis (here s and h) as 



diurnal 

yar 

seasonal 

yar 


dv 

dt 

dv 

dt 


Ry (s,cp)i, 

Ry (h, Ob 

Ry(K, O r y(s, 


( 10 ) 

( 11 ) 


( 12 ) 


where the R are the following general rotational matrices 


Ry(s, </>) = cos</>I + sin</>Ri(,s) + (1 — cos^) R2(0> (13) 

Ry(K, 0 = cosP — sin£Ri(/Z) + (1 — cos0®2(0 (14) 


with 
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o — s z 


Ri(s) = 


VSI + S 2 y + Si \ 


Sx 

0 


(15) 


S X Sy S X Sz 


Ra(l) = 


&x Sy “t - s\ 


0 X 0y Sy 

S X S Z SySz 


(16) 


Ml (h) = 


0 

h z 


\J h l +hl+hl \ _ h 


—h z 

0 

h x 


hy 

h x 

0 


(17) 


hl+hl + hi 


h x h z 


hx hy 


hyhz 


(18) 


In the specific ID heat conduction model of iBro 1 (l200fih . 
the angles have the following explicit closed forms 


tan <f> = 



( 20 ) 


where 0 < Q yar < 1 is a value satisfying Q abs = Q re f + Q yar, a is 
the Stefan-Boltzmann constant, e is the target’s emissivity, and 
the orbital period II can be expressed in terms of the target’s 
mean motion n as II = 27r/n = 2im 3 ^ 2 /VG (M + m). The an¬ 
gles <f> and £ can reach physically meaningful values of between 
0° and 45°; for more de tails on this result, and the derivation 
of equations (fl9l20l) , see IBro 1 (1200611 . 


2.6.3 Yarkovsky magnitude 


Having established the direction of the Yarkovsky acceleration, 
we now determine the corresponding amplitude. This amplitude 
may be expressed with momentum p, power P and temperature 
T using the relativistic energy momentum equation for photons 
with negligible rest mass as 



1 dp A P 
m dt me 


-±- [ P(Ti) - P(T 2 )}. 


( 21 ) 


Equation m emphasizes that the Yarkovsky effect exists 
only in the presence of a thermal imbalance on the target’s sur¬ 
face. The magnitude of the imbalance is model-dependent. We 
provide insight into our model by considering the incoming and 
outgoing power. By assuming radiative balance @ we have 


6 After reaching its equilibrium temperature the target is able to 
radiate all of the incoming energy eventually, so that it is no longer 
heating up. 


T(t) AinQyar 
4nr 2 ’ 

AouteaT*. 


( 22 ) 

(23) 


Equating the incoming and outgoing power gives 

T= / Ain L(t)Q yar \ 1/4 
\ A ou t 47recrr 2 ) 


(24) 


The target’s equilibrium temperature T depends on ratio 
of Ain/Aout. This ratio is the key to modelling the Yarkovsky 
effect. Consider first the general case of a fast rotating spheri¬ 
cal target with zero obliquity. Because the incoming radiation 
is directional and the emitted radiation is omnidirectional and 
symmetric with respect to the equator and the rotation axis, 
Ai n = ttR 2 corresponds to the target’s cross section, whereas 
4l ou t = 4-7T.R 2 . Alternatively, in the limit of slow rotators (e.g. 
a 1:1 spin orbit resonance) the target shows always the same 
hemisphere towards the radiation source (e.g. the Moon as seen 
from the Earth). In this case, only one hemisphere is constantly 
heated, and without an atmosphere or unusually high thermal 
conductivity of the body, only the illuminated hemisphere will 
be able to emit the incoming energy. Hence, A ou t = 2nR 2 and 
Ai n /A 0 ut = 1 / 2 . 

In the limit of fast rotation, all of the target’s surface ele¬ 
ments are heated up to the equilibrium temperature, and hence 
no strong thermal gradients, and no Yarkovsky effect, is ex¬ 
pected. In the limit of slow rotation, because only half of the 
target’s surface can be used as a radiator, the equilibrium tem¬ 
perature on the day side is higher by a factor of 2 1 / 4 than the 
equilibrium temperature of the entire target. In this case, 


P(7i) — P(T 2 ) = Ai£ct(2 1 /4 T) 4 - A 2 eaT 4 

- i ■ <»> 

where Ai = A 2 = A = nR 2 are the momentum-carrying cross- 
sections of the asteroid’s heated and non-heated hemispheres, 
respectively. However, the violent nature of the GB environment 
suggests that targets with all t ypes of spins a re expected, unless 
they are far from the star dVeras et alj|2014q ). Therefore, gener¬ 
ally, the temperature difference will be in between the extremes 
of 0 and 2 1 ^ i T. Consequently, 


(- fe AL(f)Q yar 
\ dt ) 4-Kmcr 2 

\ / yar 


(26) 


where 0 ^ k 1/4. The value of k is strongly linked to the 
targets’s rotation. For objects that spin very fast, i.e. E cr it —> 0, 
k = 0. Alternatively, for E —> n, k = 1/4. 

In the Solar system, the albedo of an asteroid is equal to 
(1 — Qyax). Typical asteroids have albedos between 0.1 and 0.3, 
meaning that usually Q yar is closer to unity than to zero for 
these objects. 


2.6.4 Final result 

We now summarize the findings in subsections 2.6.1-2.6.3 and 
can express the perturbative acceleration due to the Yarkovsky 
effect as 
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Figure 4. A visual summary of the relevant Q values in equation 
J29}. which represent the efficiency of absorbed radiation (Q a bs)j im¬ 
mediately reflected radiation (Q re f), and radiation which is reflected 
after a delay (Q yar ). 4> is the thermal lag angle. 


( - kAfff)Qy^Yt 
V dt ) firmer 2 

\ / yar 

where 

¥e1 y (s^) 1 y (^). 


(27) 


(28) 


2.7 Complete expression 

Finally we combine equations ©. 0 and (E3, along with k 
from equation (1261) , into the following single compact expression 


( 


dv\ 

dt J 


ra 


AL(t) 

Airmcr 2 


Qabsl + QrefI + &QyarY 




(29) 


where the I and Y represent 3x3 matrices (see equation 1281) . 
These matrices reveal the relevant regimes of motion. The first 
term is unchanging and applies for targets of any size. The third 
term is zeroed out for targets which are too small or do not 
spin. If those targets are smaller than the typical wavelength of 
incoming radiation, then only the middle term is affected. The 
Q quantities in equation (1291) are all characterized visually in 
Figure [1] 


3 RELATIVE IMPORTANCE OF DIFFERENT 
ACCELERATIONS 

Now that we have developed an expression that includes all of 
the effects from radiation, we can estimate the relative impor¬ 
tance of radiation compared to mass loss, wind drag and simply 
unperturbed Keplerian motion at a given time. First, for ease 
of notation, define the term in brackets from equation 01 as 

Q = Qabsl + Qref I + kQy ar Y. (30) 

We now describe the physical interpretation of this model. First 
the asteroid absorbs the momentum of the incoming light with 


an efficiency of Q a bs- Then part of the radiation is either re¬ 
flected immediately with an efficiency Qref and/or thermal- 
ized. In the latter case the factor kQ yal determines the relative 
strength of the Yarkovsky drift. 

Because for large bodies (A <C d), 0 ^ Q&hs ^ 1 and 
0 ^ Qref ^ 1, each component of Q must lie between 0 and 
2 inclusive. Generally, Q ^ 2 corresponds to large particles like 
asteroids only as smaller particles can retain scattering (and 
thus extinction) cross sections much larger than their geometric 
cross sections. The Yarkovsky effect has the most importance 
when Q re f = 0. In this case only a quarter of the total remitted 
radiative recoil can contribute to the Yarkovsky perturbative 
acceleration, as k St 1/4. However, this rudimentary analysis is 
deceiving, as we will show later with the averaged equations of 
motion in orbital elements. 


3.1 Radiation versus gravity 


More broadly, we can develop a sense of when radiative effects 
are important by comparing the terms on the RHS of equation 
0 - The first term is the unperturbed two-body term, the sec¬ 
ond is the acceleration due to mass loss, and the third is the 
acceleration due to radiation. Define fH as the ratio of the ra¬ 
diative acceleration to the total gravitational acceleration, i.e. 
the ratio of the third term to the sum of the first and sec¬ 
ond. The sum of the first and second terms is equal to just 
— G\M(t) + mlr/r 3 , i.e., the first term with a time dependence 
(IOmarovlll962l : lHadiidemetrioul Il963l ; iDerbtTt] 1 1 9831) . Then, be¬ 
cause |<7[ ~ 1 and |Q| ^ 2, 


AL(t) _ 3 L(t) 

2nGm (M ( t ) + m) c 8nGcRpM (t) ’ 


(31) 


an expression which is independent of the semimajor axis a and 
the eccentricity e. 

We plot /3 in Figure [a] for targets with radii of 1 m (top 
panel) and 1 mm (bottom panel), corresponding to large and 
small pebbles. The plots are a function of time, for the entire 
lifetime of stars with ZAMS (zero-age main sequence) masses of 
1.0, 2.0, 3.0, 4.0 and 5.0 Mg. The variations in each curve are 
due to the time dependence of both M and L, particularly on the 
GB p hases. We compute stellar evolutionary tracks from the SSE 
code dHurlev et al.ll200(il ~). assuming for all stars Solar metalicity, 
a Reimers mass loss coefficient of 0.5, and a superwind on the 
Asym ptotic Giant Branch (AGB) phase (IVassiliadis fe Woodl 
Il993ll . 

The earlier, minor peak in the curves corresponds to the tip 
of the Red Giant Branch (RGB) phase, and the later peak the tip 
of the AGB phase. The figure demonstrates clearly that 1mm- 
sized targets which remain gravity-dominated for the entire 
MS and RGB phases will suddenly become radiation-dominated 
during the AGB phase. This behaviour remains true regardless 
of the location of the targets. For larger targets which are about 
one meter in size, the contribution of radiation to their orbital 
motion along the AGB may reach the one per cent level, poten¬ 
tially enough for example to perturb the targets out of mean 
motion resona nce with ano ther object - a subject for a future 
study (see e.g. |Pastorll2014l 'l. 

These results must be considered in the appropriate con¬ 
text; they apply only for snapshots of system evolution. In Sec¬ 
tion 4.3.4 we will show that radiation may in fact dominate the 
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I | 


_ l 


K AL(t)M(t) ( dM{t) \ 1 _ 3 L 

~ 7r mcr 2 v \ dt ) 4'KRpncr 2 v'h{t) 

3 L(t) (1 + ecos /) 2 

4-n-GM{t)Rpc'li(t) (1 - e 2 ) 3/2 (1 + e 2 + 2ecos /) 1/2 

(33) 


where 


*(*) 


1 

nM{t) 


dM(t ) 

dt 


(34) 


In equation (1331) . we have introduced the true anomaly / 
and expressed the distance and velocity in terms of orbital ele¬ 
ments as 


10 * 

100 

1 

es. 


10 * 


10 4 

10 100 500 1000 5000 10* 

Time since star was born (Myr) 

Figure 5. Approximate relative importance of the acceleration due 
to radiation versus that from gravity (/3), as a function of stellar age, 
for two different pebble-like targets with radii R. The target density 
of 2 g/cm 3 corresponds to a typical density of an icy Kuiper belt 
body. Both the stellar mass and luminosity affect 0, and vary strongly 
with time during the GB phases of evolution, causing the peaks in 
the curves. Particles smaller than about 1 mm will be affected more 
strongly by radiation than gravity at some point in their evolution. 
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1 


secular evolution of targets which are orders of magnitude larger 
than one meter. The reason is due to a potentially accumulating 
drift from the Yarkovsky effect, which cannot be deduced from 
the ratio of instantaneous accelerations. 


3.2 Radiation versus mass loss 

The mass lost by the star, particularly during the GB phases, 
will contribute to the orbital motion through a change of grav¬ 
itational potential in the system. We can compare the orbital 
changes due to radiation to those due specifically to mass loss 
from equation m through 

_ 1 1 mm m 

\dt / ml 2 (M(t) + m) dt ' y 1 

Equation (El helpfully identifies the acceleration due to 
isotr opic mass loss from the star in terms of positions an d veloc¬ 
ities (lOmarovll 19621 ; lHadiidemetrioull 19631 : iDepritll 19831 1. There¬ 
fore let us define 7 as the magnitude of the ratio of the acceler¬ 
ation due to radiation to the acceleration due to mass loss. We 
obtain 


all — e 2 ) 

r = — -— 

1 + e cos / ’ 

v = , \/l + e 2 + 2 ecos f. 

VT^e 2 


(35) 

(36) 


The dimen s ionles s mass loss adiabaticity index 'k is from 
IVeras et al.l (l201ll l. This index provides a scaled ratio of the 
target’s orbital period to mass loss timescale, and distinguishes 
two distinct regimes of motion: adiabatic (4/ -C 1), when 
the semimajor axis increases at a constant rate, and run¬ 
away ('k 1 ), where o, e, and the argument of pericen- 

tre w, increase or decrease at nonuniform rates. This index 
can be used to det ermine the evolution of individual known 
planetary systems dVeras fc Wvattl 120121 : iMustill et al.1 120131 : 
iTadeu dos Santos et al.ll2015l) , isolat e the strength of thre e-body 
interactions amidst stellar mass loss dVovatzis et alj|2013f ). eval¬ 
uate the contribu tion of escaped planets to the free-floating 
planet population dVeras fc Ravmondll2012| l. and determine the 
relative i mportance of other perturbative effects such as Galac¬ 
tic tides (jVeras et al.112014(11 ). In equation (1331) and throughout 
the paper, we assume that the mass loss is isotropic, an excellent 
approx imation for targets within a few hundred au dVeras et al.1 
l2013bll . 

Adopting limiting values of / help us understand equation 
El- When the target is at pericentre, the term in square brack¬ 
ets becomes (1 + e)” 1,/2 (l — e)“ 3,/2 ; at apocentre, this term is 
(1 — e)~ 1 ,/ 2 (l-|-e)~ 3,/2 . Therefore, for highly eccentric targets, the 
acceleration from radiation becomes much more important than 
that from mass loss, more when the target is at pericentre than 
at apocentre. For moderate or small eccentricities, the brack¬ 
eted term may be approximated as unity, as for e = 0.5, the 
term approximately equals a value between 0.8 and 2.3. Also, 
unlike /?, 7 is dependent on the semimajor axis, through 'k(f). 
More generally, the importance of radiation relative to mass loss 
increases for (i) adiabatic motion, (ii) smaller targets, (iii) less 
dense targets, and (iv) less massive stars. 

We chart 7 in Figure El for small exobody analogues of the 
asteroid belt (left panel) and Kuiper belt (right panel). The 
left and right j/-axes in both plots correspond to metre-sized 
and millimetre-sized targe ts. We compute the mass loss index 
in the same manner as in IVeras fe Tout! J 2012 f) and assume ^ 
is defined in terms of n rather than II. This difference has a 
negligible impact on the resulting curves in the figure, which 
span several orders of magnitude. We consider the targets at 
their pericentres; if instead we considered their apocentres, then 
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Radiation / mass loss 
for exo-Asteroid belt-1 i ke target 



Radiation / mass loss 
for exo-Kuiper belt-like target 
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Figure 6. Relative importance of the acceleration due to radiation versus mass loss ( 7 ), as a function of time, for a metre-sized target (left 
y-axes) and a mm-sized target (right y- axes) in an exo-Asteroid belt (left panel) and in an exo-Kuiper belt (right panel). The variables e, a, / and 
p refer to the target’s eccentricity, semimajor axis, true anomaly and mass density, respectively. In all cases, orbital evolution due to radiation is 
more important than that from stellar mass loss, even at the tip of the AGB phase. 


the curves would be shifted downward by factors of 2/3 (left 
panel) and 3/7 (right panel). 

The curves suggest that for metre-sized and smaller targets 
within tens of au, radiation is much more important than mass 
loss for orbital evolution. Although this ratio is most extreme 
on the MS, neither radiation nor mass loss contribute much to 
the unperturbed Keplerian motion (given by the first term on 
the RHS of equation [lj during this phase. On the AGB phase, 
this ratio is minimized, but radiation is still more important 
than mass loss. Therefore, when modelling sub-asteroidal-sized 
targets along the GB phases, we assert that considering only 
their orbital evolution due to mass loss is insufficient. 


the gas molecules. Recent applications which partly rely on re¬ 
sults from physical experiments do not all use the same val¬ 
ues f or the different reg imes (see equations 2 -4 of lGaraud et al.l 
12004, equations 15-20 of Gibbons et all2012l . and equations 6-11 
of Loren-Aguilar fe Bata 20141) . demonstrating slight disagree¬ 
ment. _ 

We use the equations from lGaraud et ahl d2004l ). which can 
be expressed in our notation as 



te) (%-*)> ^«c 

(7^) (v e -v)\v e -v\, R »C 


(37) 


3.3 Radiation versus wind drag 

In the previous subsection, we focused on the gravitational evo¬ 
lution due to mass being lost from the two-body orbit. We con¬ 
sidered the mass lost to be through a gaseous stellar wind, and 
ignored the time lag from when the star loses mass to when the 
ejecta reaches the target. Such subtleties may be unimportant 
for targets the size of planets, but can be crucial for smaller 
targets. 


where p g is the density of the gas, ( is the mean free path length 
of the gas, v g = (v gx , v gy , v gz ) is the velocity of the gas, v s is the 
local sound speed, and 



0.165, 


Re ^ 1 

1 < Re < 800 ( 38 ) 

Re > 800 


3.3.1 The different regimes of wind drag 

A solid moving through a gaseous medium will experience a 
drag force. This well-studied force, in its most basic form, is 
proportional to the square of the solid’s velocity. This form has 
imp ortant applications f or protoplanetary discs (e.g. equation 
3 of llwasaki et al.l 2002|, equations 3-4 o f Beauge et all l20ld . 


and equation 1 of lYamaguchi fe Kimural 120141 ) but is applica- 

ble to a diverse range of systems, including for example the 
the current P luto-Cha ron system environment (equation 3 of 
[Pprter_^Grundv 2014) and the interstellar medium (equation 
3 of Howe fe R afikov 20 li|). However, the drag coefficient is, in 
general, not constant. Equations 3.5-3.10 of lAdachi et al.l Jl976l ) 
show how the coefficient is actually a function of the Reynolds 
number, the Mach number and the mean free path length of 


the Reynolds Number Re can be expressed as 
\v s ~ v\. (39) 

The two pieces of equation G3> refer to the Epstein regime 
(upper relation) and the Stokes regime (lower relation). 

Equations m-m should be applicable to GB star winds, 
just as they would be f or protoplanetary discs. In their study 
of drag in GB systems, iDong et alj d20ld) utilize one form of 
the drag equations in the Stokes regime, and for their model 
suggest that the drag contribution is much smaller than grav- 
ity in the radia l, but not necessarily tangential, direction. 
iBonsor fe Wvattl (l2010i ) instead compute a spiral-in timescale 
due to GB stellar wind drag. 


such that 


Re=^ 
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12 Veras, Eggl, & Gansicke 

3.3.2 The velocity lag 

A key quantity in equations (I37I) - (I39D is the difference in velocity 
between the gas and target. This difference is also proportional 
to the Reynolds number (equation [39]). Therefore, let us explore 
this gradient in more detail. Each component of both v and v g 
can be expressed in terms of the target’s orbital elements. For 
v, we use equations 5-6 of I Veras fe Evansl ]2013al ') and define a 
rotation matrix R with components 


R 11 

= cos fl cos u> — sin Q sin lo cos i, 

(40) 

R 21 

= sin S2 cos lo + cos fl sin a; cos i, 

(41) 

R 31 

= sin lo sin i, 

(42) 

Rl2 

= — cos 12 sin lo — sin cos lo cos i, 

(43) 

R 22 

= — sin fl sin lo + cos cos lo cos i, 

(44) 

R 32 

= cos lo sin i, 

(45) 

Rl3 

= sinflsini, 

(46) 

R 23 

= — cos Id sin i, 

(47) 

R 33 

= cos i, 

(48) 


where i is the inclination of target’s orbit, and and lo are its 
longitude of ascending node and argument of pericentre. The 
target’s velocity components are then 


/ 

v = R 

V 


— na sin / 


na(e -\-cos /) 



o 


\ 

) 


(49) 


The velocity of the gas (or wind) is assumed to be constant 
and in the radial direction. This velocity is dictated entirely 
by the star, but shares the same space, and geometry, of the 
target. Consequently, we may express the components of v g at 
the target’s location in terms of the orbital parameters of the 


target as (equation 2.122 

of Murrav & Dermott 1999li 


U gx = \V B \ 

[cos n cos (lo 

+ /) 

— sin S2 sin (lo + /) cos i] , 






(50) 

v sy = 1 1 

[sin 12 cos (lo 

+ /)■ 

-1- cos Q. sin (lo + f ) cos i ], 






(51) 

U gz = |Tg| 

[sin (lo + /) 

sini] . 


(52) 


Equations m-(59 state that the gas velocity is in the radial 
direction. However, we must carefully consider the context when 
applying equations ([50])-(52]). For example, they would be incon¬ 
gruous if inserted into averaged equations of motion. 

The magnitude of v g is dictated by the physical properties 
of the sta r. Because the launch ing mechanism of the wind is 
unknown (Bladh fc Hofnei 2012)) a nd may vary as a function of 
stellar evolution isee lOwoc ki l2004l for a review), we can make 
only an approximation. We appr oximate the e jecta speed as a 
function of the star’s latitude as llOwoc ki 120131") 


(2GM\ 


1 - 


GM 


cos 2 Z 


(53) 


where f?* is the radius of the star, Z is latitude, and Kot is the 
rotational velocity of the star at the equator. Time dependencies 


are included to emphasize the potentially significant variance of 
these quantities along the GB phases. MS stars do also have 
winds, but for fixed values of f?* and Vrot- WD stars have not 
yet been observed to emit winds, and therefore in this context 
the concept of gas drag through stellar winds does not exist and 
equation (53]) is not applicable. 


3.3.3 The gas density, mean free path length, and sound speed 


Before proceeding, we must make some approximations about 
the gas density p g , the mean free path length (which we denote 
as £), and the sound speed v s . If the wind is spherically sym¬ 
metric so that p g can be deter mined from the m ass loss rate, 
then we can write (equation 5 of lDong et al.l[201(i| j 


dM(t) 

dt 


47 rr 2 p g |u g | . 


(54) 


With a given mass loss rate, and a gas velocity given by equation 
(1551) . one can determine p g . 

Given an expression for p g , how can we determine (j? These 
quantities are related through 


c 


pm-H 

OPs ’ 


(55) 


where p is the mean molecular mass of the gas (or stellar wind), 
niH is the mass of a Hydrogen atom, and 6 is the collision cross- 
section of the gas molecules. In equation (1551) . rnn is fixed and 
p is on the order of unity, so we need to bound the collisional 
cross-section 9. Based on the known radii of chemical elements, 
10“ 2 O m 2 < 9 < 10~ 18 m 2 . Consequently, we make the approxi¬ 
mate relation 


p g C~10 8 kg/m 2 . (56) 

By using this relation, we plot £ in Figure0 The plot illustrates 
how at a fixed point in space, £ will vary by many orders of 
magnitude during the GB phase. £ is minimized at the tip of 
the AGB, and scales with the square of the distance to the star. 

The sound speed v s is also u nknown, and is a function of the 
properties of the changing wind. lOwockil 420141 1 suggests that v s 
is a small fraction (~ 0 . 001 ) of the escape velocity of the wind 
from the star. For simplicity and for lack of better constraints, 
we will assume that v s is constant and given by this ratio. 


3.3.4 The Reynolds number 

Having now established methods or values for obtaining <4 and 
v s , we can attempt to estimate the Reynolds number (equation 
1391) . The Reynolds number is particularly important because 
it defines the flow regime in which the target resides. Figure 
[H] demonstrates how the variation in Re can be great enough 
for a fixed point in space to be subject to several different flow 
regimes throughout GB evolution. 


3.3.5 Comparing wind drag with radiation 

Having obtained some physical intuition for the variation of the 
Reynolds number, we can now finally consider the acceleration 
due to wind drag. In the same vein as we have defined (5 and 7 , 
here let us define 5 as the ratio of the perturbative acceleration 
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Radiative and drag forces in post-MS systems 13 


Mean Free Path Length Evolution 



Figure 7. Evolution of the mean free path length of the stellar wind 
for a fixed-in-space target across the giant branch phases of stellar 
evolution. The target’s physical and orbital parameters are R = 1 
cm, a = 5 au (left y-axis), a = 50 au (right y-axis), e = 0.2, i = 0°, 
= lj = f = 0°. The plot demonstrates how £ undergoes changes of 
over 5 orders of magnitude, and achieves a minimum at the tip of the 
AGB. 


Radiation/wind drag 



Figure 9. Relative importance of radiation versus wind drag for a 
fixed-in-space target across the giant branch phases of stellar evolu¬ 
tion. The R = 1 cm target has orbital parameters of a = 5 au, e = 0.2, 
2 = 0°, and = lj = f = 0°. However, the evolution of S is largely in¬ 
sensitive to any of these parameters and is instead largely determined 
by the competition between the mass and luminosity evolution of the 
star (equations |57| and l58b . 


GB Reynolds Number Evolution 
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Figure 8. Evolution of the Reynolds number of the stellar wind for 
a fixed-in-space target across the giant branch phases of stellar evo¬ 
lution. The target’s orbital parameters are a = 5 au, e = 0.2, i = 0°, 
Q = cj = f = 0°. The left and right y-axes respectively illustrate the 
evolution for a R = 1 cm pebble-like target and a R = 1 km asteroid¬ 
like target. The asteroid traverses all three of the Stokes flow regimes 
(equation l38(l . 


due to radiation to the perturbative acceleration due to wind 
drag. In the Epstein regime, 


3 L (1 + ecos /) 2 
8nca 2 (1 — e 2 ) 2 p g |u g | 


3 L ( dM \ -1 

2 c |{?g — v| \ dt J 

whereas in the Stokes regime 


(57) 


5 < 


3 L (1 + ecos /) 2 ._ 
87 rca 2 (1 — e 2 ) 2 pg-B s 
3L |v g | / dM \ 
2cB\v e -v\ 2 \~df J 


(58) 


The simplification we performed in equations (|571l - (|58l> 
holds only when one assumes the spherically symmetric wind 
of equation (1541) . In this approximation, 5 is largely indepen¬ 
dent of the target’s location or size, and instead is dependent 
almost entirely on the properties of the star. Figure [9] illustrates 
the evolution of 5 for the stellar models that we adopted. The 
curves suggest that the competition between radiation and wind 
drag is close, and both should be considered together in GB 
planetary studies. 


3.4 Comparison of all forces together 

Having defined the dimensionless ratios /?, 7 and 5, we can now 
compare the strength of mass loss, gravity, wind drag and ra¬ 
diation at given snapshots of time. See Figure [10] In all cases, 
radiation and wind drag are comparable in strength, and com¬ 
parable at all target radii sampled. Recall that S is independent 
of R everywhere except when the target is in the Stokes regime 
and Re < 800. However, the pebble (with Re ~ 10 ~ 9 —10~ 2 ) and 
asteroid (with Re ~ 10 -4 — 10 3 ) are both always in the Epstein 
regime here, and when the planet (with Re ~ 10° —10 7 ) is in the 
Stokes regime, usually Re > 800. Consequently, changes in the 
purple dotted fines are almost indiscernible. Expectedly, gravity 
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dominates the evolution of planet-sized targets, and mass loss 
becomes relatively more important the further one moves away 
from the star. These trends vary little for different progenitor 
stellar masses. 

Although broadly useful, this type of comparison does not 
(i) characterize the evolution of individual orbital elements (such 
as eccentricity), and (ii) fails to account for cumulative effects. 
Consequently, the middle panels suggest that gravity will dom¬ 
inate the evolution of asteroids. These plots do not reveal that 
the long-term effect of the (radiative) Yarkovsky force can ex¬ 
cite an asteroid’s orbit to instability, as shown in Section 4.3.4. 
These plots are primarily useful for order of magnitude estimates 
and short-term evolution; detailed evolutionary models should 
instead rely on Figure [T] 


4 EQUATIONS OF MOTION IN ORBITAL 

ELEMENTS 

The physical intuition gained in the last section from the com¬ 
parison of forces provides a foundation for the more precise 
treatment that we supply here. We seek to express the equa¬ 
tions of motions due to wind drag and radiation entirely in 
terms of orbital elements t hrough the formal ism first men¬ 
tioned in IVeras et alj (1201 ill and developed in IVeras&Evansj 
(|2013a|). Based on work bv lEfroimskvl (l2005blf and lGurhll ( 20071 4 . 
IVeras fe Evans! (l2013al l delineated how one can derive orbital- 
element-only equations from a two-body problem with a pertur¬ 
bation expressed entirely in terms of Cartesian elements, with¬ 
out any averaging. This procedure assumes only boundedness of 
the orbits, and otherwise does not make any assumptions about 
the magnitude of the perturbation or the resulting orbital ele¬ 
ments. 

Among the benefits of this task is the ability to then average 
over the result to create a new set of equations which can be 
integrated much more quickly and provide perspective on the 
target’s secular evolution. Further, for the particular case of 
perturbations due to radiation, we will demonstrate that in some 
cases averaged leading-order terms vanish, helping to provide a 
better assessment of the orbital changes than those from Figs. 
[5] and [6] Henceforth, for ease of notation, we drop the explicit 
dependence of time on M and L. 

Unlike the other variables, the time evolution of the true 
anomaly need not be split up by type of force. The following 
relation holds true individually for mass loss, wind drag and 
radiative forces, and together for the overall evolution 


df _ n (1 + ecos f) 2 
dt ~ (1 — e 2 ) 3/2 


duj 


. dUl 


, — cosi—— 
dt dt 


(59) 


The true anomaly also represents the variable over which we 
average in order to obtain the secular equations. 


4-1.1 Unaveraged equations 

The unaveraged equatio ns of motion i n orbital elements for 
isotropic mass lo ss are l|Omaro 3 ll962j ; lHadiidemetrioul [l963i ; 
IVeras et al.ll201lh 



a (l + e 2 + 2ecos/) 1 dM 
1 — e 2 M + m dt ’ 


- (e + cos /) 


1 _ dM 

M + m dt ’ 


(60) 

(61) 


0 , 


(62) 


0 , 


(63) 


sin / 1 dM 

e M + m dt 


(64) 


The m ore complex anisotropic equations of motion (IVeras et al.l 
l2013lJl in fact more realistically represent the motion, but, as 
previously suggested, they provide a negligible improvement on 
equations (I60II64I) except in extreme cases. 


4-1-2 Averaged equations 


We obtain averages by performing the following integral for each 
variable (here for averaged semimajor axis change due to mass 
loss) 



TfC-) df- 

2n Jo \ rft /mi (1 + ecos/) 2 


(65) 


We find that all secular motions vanish except for the semimajor 
axis’ 


/ / da\ \ _ a dM 
\\dt) ml / M + m dt 



( 66 ) 

(67) 


The averaging assumes that dM/dt is a constant and is much 
less than df /dt, such that M(t) is also approximately constant 
on orbital timescales. Equations l[66|l-(|67|l represent an excellent 
appr oximation for targets within hundreds of au (IVeras et al.i 

ml ])• 


4.2 Wind drag equations 

Our derivation of the wind drag equations of motion in orbital 
elements make no assumption about v s = (v sx ,v gy ,v gz ). The 
gas velocity here does not have to adhere to the prescriptions in 
equations (I50l) - (l53l) nor equation d54l) . 


4.1 Mass loss equations 

First, for completeness and comparison, we repeat the known 
isotropic mass loss equations. Mass loss significantly changes 
the orbits of targets of all sizes. 


4-2.1 Auxiliary expressions 

In order to express our equations in a compact manner, we define 
a standard set of auxiliary variables that appear regul arly in 
studies of the perturb ed two-body pr oblem (IVeras et al.ll2013bl : 
IVeras fe Evandl2013al lbl: IVeradl2014al l . 
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Comparison of all effects 
(dot-dashed = 1//3, dashed = 1 / 7 , dotted = 1/5) 
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Figure 10. Relative instantaneous acceleration strengths of gravity (1//3; orange dot-dashed lines), mass loss (I/ 7 ; blue dashed lines), and wind 
drag (l/<5; purple dotted lines), all compared to radiation for Mzams — 2 Mq stars along the GB phases. The relative strengths and profile 
features for other progenitor stellar masses are similar. The top, middle and bottom panels illustrate targets which are a planet, asteroid and 
pebble, and the left and right panels represent locations given by a = 5 au and a = 500 au assuming the target is on a circular coplanar orbit. 
This plot represents just a snapshot in time and does not take into account potentially large long-term accumulations; detailed models should 
instead consider all checked forces from Tabled 
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Ci = e cos uj + cos (/ + ui) 

C 2 = e sin uj + sin (/ + u) 

C3 = cosi sin £2 sin (/ + uj) — cos £2 cos (/ + uj) 
C4 = cosi cos£2 sin (/+ cj) + sin£2cos (/+ w) 
Cs = (3 + 4e cos / + cos 2/) sin uj 

+ 2 (e + cos /) cos uj sin / 

Ce = (3 + 4e cos / + cos 2 /) cos uj 

— 2 (e + cos /) sin uj sin / 

C 7 = (3 + 2e cos / — cos 2/) cos w + sin ui sin 2/ 

Cs = (3 — cos 2/) sin a; — 2 (e + cos /) cos w sin / 

Cg = (3 + 2ecos / — cos 2/) sin w — cos wsin2/. 


( 68 ) 


Additionally, in the following derivations we found that a 
commonly-occurring quantity is 


n = 


4 “ ( + 


U S x + 


an (C 2 cos £2 + Ci sin £2 cos i) 


an (—Ci cos i cos fi + C 2 sin £2) 




, anCi sin i 

+ 1 Wgz ~ 


1/5 


(69) 


4-2.2 Unaveraged equations for Epstein regime 

Recall that in the following derivations, we make no assumptions 
about v g . Consequently, from equation cm we obtain 


+ 


2u a pg 


nRp (1 — e 2 ) 


— an (l + e 2 + 2e cos /) 


Vi -> 


VgzCi sin i + Vg y (C 1 cos i cos Q. — C 2 sin Q) 


— u gx (Ci cos i sin Q + C 2 cos f2) 


(70) 


de\ 

dt J , 2anRp (1 + e cos /) 


VsPg 


4aen 


+4an cos / (l + e 2 + e cos /) + Vl~e 2 
+v gy (—Ce cos i cos Q. + C5 sin S2) 

+ v gx (Ce cos i sin + C 5 cos Q.) 


— Vg Z Ce sini 


(71) 


Vi — e 2 v s pg cos (/ + uj) 


di 


dt J dr anRp (1 + e cos /) 


x \vgz cos i — Vgy sin i cos Q + v gx sin i sin $7] . 


(72) 


dQ, 

dt 


Vi — e 2 v s p g sin (/ + uj) 
anRp (1 + ecos /) 


x [v gz cot i — Vgy cos fi + Vg X sin fi], 


( 73 ) 


VsPg 


AanRpeVi — e 2 (1 + e cos /) 


8 anV 1 — e 2 (1 + ecos /) sin / — 2 (l — e 2 ) 
— Vgz [Cg sin i + 2e cos i cot i sin (/ + w)] 
—Vgy (Cs cos i cos O + C 7 sin fi) 

+ Vgx (— C 7 cos O + Cs cos i sin S 2 ) 


(74) 


Equations G2D-GD reveal that all the orbital elements vary 
along a single orbit due to Epstein drag, even if the orbit is ini¬ 
tially eccentric or circular. This fact remains true even if spher¬ 
ical symmetry is imposed. 

4-2.3 Unaveraged equations for Stokes regime, Re 4 1 

The orbital element evolution equations of motion in this regime 
are equivalent to those in Section 4.2.2, but with right-hand- 
sides multiplied by a factor of 2iC,/2R. Recall that targets will 
rarely find themselves in this regime, and must be close to the 
star to do so. 


4-2-4 Unaveraged equations for Stokes regime, 

IV Re V 800 

The equations of motion in this subregime have a similar form 
to those in Section 4.2.2, but are scaled differently. 

da\ _ 3 ■ 6 2/5 C 3/5 u a 3/5 p E ^ 

dt J dr Vanp(i-e2)R*/ 5 


: jna 3 / 2 (l + e 2 + 2ecos /) + aV 1 — e 2 


— VgzCi sin i + Vgy (—Ci cos i cos Q + C2 sin fi) 


+ Vgx (Ci cos i sin fi + C2 cos fi) >, (75) 


de\ _ 3 7 / 5 C 3 / 5 u a / 5 p g -fc 

dt J di 2 8 / 5 a 3 / 2 np (1 + ecos /) R s / 5 

x 12a 3 ^ 2 n [2 cos / (l + e 2 ) + e (3 + cos 2/)] + aV 1 — e 2 


— VgzCe sin i + v gy (—Ce cos i cos Q + Cs sin £2) 


+ Vgx (Ce cos i sin £2 + C5 cos £2) 


(76) 
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di\ _ 3 7 ' /5 C 3 ^ 5 fs' /5 p g 7?.Vl — e 2 cos (/ + uj) 


dt) 


2 3 / 5 anp (1 + e cos /) i ? 8 / 5 


x < Vgz cos i — Ve Y sin i cos £2 + Vex sin i sin £2 


d£ 2 \ _ 3 7 / 5 £ 3 / 5 Ug/ J p s TZy/l — e 2 sin (/ + oj) 

dt ) dr 2 3 / 5 anp (1 + e cos /) R s / 5 

x ^ v S z cot i — v sy cos £2 + v g x sin SI j 


dio 


3 7 / 5 C 3/ Se 3/5 p E ^ 


dt ) dr 2 8 / 5 anes/\ — e 2 p (1 + ecos /) R 8 / 5 


4an\J 1 — e 2 sin / (1 + e cos /) — (l — e 2 ) 

— Vgz (Cg sin i + 2e cos i cot i sin (/ + a;)) 
—Vgy (Cs cos i cos £2 + CV sin £2) 

+ Vgx (Cs cos i sin £2 — C 7 cos £2) 


(77) 


(78) 


(79) 


di\ _ K.pgVJ’^y/l — e 2 cos (/ + a;) 
dt J dr anRp(l + ecos/) 


x v gz cos i — Vgy sin i cos £2 + v gx sin i sin £2 }, 


dn\ _ K.pg'R 5 l 2 VT^e I sin (/ + w) 
dt 


anRp (1 + ecos /) 


x ^ v gz cot i — Vgy cos £2 + v gx sin £2 } , 


du 


K,pg n 5/2 


dt ) dr 2naRpe\/\ — e 2 (1 + e cos /) 


(82) 


(83) 


- an\J 1 - e 2 (CzCt - C 1 C 9 ) - (l - e 2 ) 


n gz (C 9 sin i + 2e cos i cot i sin (/ + w)) 


—Vgy [— cos * cos £2 (Cg — 2e sin (/ + uj)) — C 7 sin £2] 


+ Vgx [— cos i sin £2 (C 9 — 2e sin (/ + uj )) + C 7 cos £2] 


. (84) 


4-2.5 Unaveraged equations for Stokes regime, 800 ^ Re 

In this subregime, B is a constant. Let B = 0.165 = 1C. Although 
the equations simplify and (di/dt) a r and (dCt/dt)di vanish, the 
remaining equations are still not generally solvable analytically. 
The full equations are 


2tCp g lZ 5/2 


da 


dt J dr nRp (1 — e 2 ) 


an (l + e 2 + 2 ecos /) 


+Vl^> 


— Vg Z Ci sin i + Vgy (—Ci cos £2 cos i + C 2 sin £2) 


+ Vgx (Ci sin £2 cos i + C 2 cos £2) 


(80) 


de\ 


K. P gR?l 2 


dt) dr 2naRp (1 + e cos /) 
x < 2 an [3e + 2 (l + e 2 ) cos f + e cos (2/)] 


+ \J 1 — e 2 — Vgz Cg sin i 
+u gy (—C 6 cos £2 cos i + C 5 sin £2) 
+ Vgx (Cs sin £2 cos i + Cs cos £ 2 ) 


(81) 


4-2.6 Averaged equations 

The unaveraged equations demonstrated exactly how the oscu¬ 
lating orbital elements change at a given point in time, and most 
practically, along a single orbit. Averaging over these equations 
would require additional assumptions, which we do not make 
here. One can, for example, assume a spherically symmetric 
wind (as in equation 1541) and then some particular analytical 
prescription for stellar mass loss to treat the variation of Vg. 
However, even with those assumptions, performing averaging 
over variables such as TZ (applicable in the Stokes regime, with 
Re > 1) becomes computationally prohibative. 


4.3 Radiation equations 

The equations of motion due to radiative effects (from equation 
1291) are too complex to be reasonably expressed in osculating or¬ 
bital elements due to the Yarkovsky contribution. Because the 
diurnal Yarkovsky acceleration contains an explicit dependence 
on s, one must also treat the evolution of the target’s spin. We 
can do so through equation m only if we assume a slow rate 
of change. Further, because the seasonal Yarkovsky contribu¬ 
tion is an explicit function of the specific angular momentum 
through equations (HI, 113, HI and (HI, the final result is 
unmanageable even if the diurnal contribution is neglected. 

Nevertheless, our perturbation method is sufficient for the 
purposes of this paper, as we aim to place bounds on the mo¬ 
tion and not model any particular system. Consequently, when 
deriving the perturbation equations below, we assume that the 
matrix elements of Q are independent of position and velocity. 
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4-3.1 Auxiliary expressions 

Even with this assumption, the final resulting equations are 
long. However, we have discovered that they can be expressed 
compactly by using a smart choice of auxiliary variables. First 
we define two additional C variables that represent commonly 
occurring quantities in the equations of motion of our problem. 


Cio = — 3esincj — 2 sin (/ + uj) + esin ( 2 / + lj), 

C n = — 3ecoscu — 2 cos (/ + uj) + ecos (2/ + u>). 

(85) 

Also, we found that one quantity which appears in all of 
the terms that are proportional to 1 /c is 

/ Qii \ / Q12 \ 

D c = — Cs I Q21 I + C4 I Q22 I 

V Q31 J \ Q32 ) 


l Q 13 \ 

+ sin i sin (/ + uj) I Q 23 ] ■ ( 86 ) 

\ Q 33 J 

A quantity appearing in all the terms that are proportional to 
1 /c 2 is 

/ Oil \ 

D c 2 = — (Cio cos 12 + C 11 cos i sin 12) Q 21 

V Q31 / 


/ Q12 \ 

— (Cio sin 12 — C 11 cos i cos 12) Q 22 

\ Q32 / 


+ Cn sin i 


(87) 


Terms specific to the orbital element expressions are 

( —Ci cos 12 — Ci sin 12 cos i \ 

—C 2 sin 12 + Ci cos 12 cos i 1, 

Ci sin i ) 


C 5 cos 12 — C 6 sin 12 cos i 
C 5 sin 12 + Ce cos 12 cos i 
Ce sin i 


D e = - 


D n = 


Di = sin iDn , 

( — C 7 cos 12 + Cs sin 12 cos i 
—C 7 sin 12 — Cs cos 12 cos i 
—Cg esc i + Cs cos i cot i 



( 88 ) 

(89) 

(90) 

(91) 

(92) 


4-3.2 Unaveraged equations 

Equipped with the above auxiliary variables, the complete final 
equations of motion entirely in orbital elements are 

fda\ _( 1\ AL (1 + ecos /) 2 - 

\dt) n \ c / 27rmna 2 (1 — e 2 ) 5/2 \ c L> ' 


+ 



/ de 
\ dt 


ra 



AL (1 + ecos/) 
47 rma (1 — e 2 ) 3 


(p c 2 • Da^j 


A ALg + ecosf) 
c J 871 mna 3 (1 — e 2 ) ' ' 

AL (1 + ecos /) 

167ima 2 (1 — e 2 ) V J 


(93) 


(94) 


/ di 
\dt 


ra 






1 \ AL (1 + e cos /) cos (/ + uj) 


(-D c • A) 


471 mna 3 (1 — e 2 ) 3 ^ 2 
AL (1 + ecos/) cos {f + uj) tg ^ g \ 
8nma 2 (1 — e 2 ) 2 ' c > 


1 \ AL (1 + e cos /) sin (/ + uj) 


(d c • Dn) 


4irmna 3 (1 — e 2 ) 3 ^ 2 
AL(l + ecos /) sin (/ + uj) (g 2 g \ 
871 ma 2 (1 — e 2 ) 2 v c > 


AL (1 + ecos/) 


c J 8iumna 3 e (1 — e 2 ) 3 ^ 2 
AL (1 + ecos /) 


167ima 2 e (1 — e 2 ) 


{p 0 ■ £>p 

2 p c 2 ■ £>P . 


(95) 


(96) 


(97) 


Equations (I93I) - (I97I) are convenient because all of the D auxiliary 
variables, which are on the order of unity, are isolated. They also 
show that the amplitudes of the orbital element rates of changes 
are independent of i, and 12. Further, for many purposes, the 
1 /c 2 term may be neglected. 


4-3.3 Unaveraged equations with no Yarkovsky 


When the Yarkovsky effect is “off”, which typically occurs for 
targets with diameters less than 1 cm -1 m, a significant simplifi¬ 
cation occurs. Now, Q is a diagonal matrix, and we can convert 
that into a scalar. Consequently, equation (12911 becomes 



AL (Qaba + Qref) ^ 

471 mcr 2 


(98) 


Subsequently, the equations of motion simplify to the following 
relations 



AL (Qaba + Qref ) (1 + e COS /) 2 
271772 (1 — e 2 ) 3 


eVl — e 2 sin / /1 


—2 — 3e 2 — 4e cos / + e 2 cos (2/) / 1 \ 

+ 2a VcV 

AL(Qaba + Qref) (1 + e COS f) 2 
471772 (1 — e 2 ) 3 / 2 



(99) 
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x 


sin / ^ 1 ^ ^ —4 cos f + e (cos (2/) — 5) / 1 
na 3 \c/ 2a 2 \/l - e 2 \ c2 


( 100 ) 


/ di 
[di, 


) 


ra 


Y=0 




AL (Qabs + Qref) (1 + e COS /) 2 
Annie (1 — e 2 ) 3 ^ 2 


( 101 ) 


COS/ f 1^ + sin/(2 - ecos/) ( 
na 3 \c/ a 2 y/l - e 2 \c 2 J 


( 102 ) 


For pebbles which are subject to intense radiation, equa¬ 
tions 11981) - (11021) describe their evolution due to this perturba¬ 
tion exactly. For asteroids, however, equations (I93M97I I should 
be used instead. For particular models which require a detailed 
analysis of the perturbation at the pericentre or the stationary 
points of the motion along an orbit, equations ®-GS2D are 
compact enou gh to allow for an analysis similar to the one by 
IVerad (l2014bll . 


4-3.4 Averaged equations 

Now we return to the general equations in (|93l)-(|97|l and average 
over them. Integrating these equations yields long expressions. 
We write out these expressions in the appendix (with the 1/c 
term only) because they can be cast in a revealing form. 

Just by inspection of equations (lA2l)-(lA6t. one can see that 
when the Yarkovsky effect does not operate (e.g. and the off- 
diagonal terms of Q become zero) all of the (1/c) are zeroed 
out. Therefore, radiative treatments which do not include the 
Yarkovsky effect are potentially missing an important perturba¬ 
tion in the system! 

Now we quantify this statement. Because the matrices in 
Appendix A are all of order unity, we can write 



1 AL \ 
c An mna 2 J ’ 

1 AL \ 
c 8nmna 3 J 

1 AL \ 
c 8nmna 3 J ’ 

a AL \ 

K c 8nmna 3 J ’ 

'1 AL \ 

. c 32nm.na 3 ) 
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(104) 

(105) 

(106) 
(107) 


Isolating the rate of change of eccentricity with time yields 


O - 


1 AL \ 0.08 




c 8 nmna 3 ) Myr \ 1 Mq ) \ 2 g/cm 


- 1/2 


( -^V 1 (—Y 3/2 (—)■ 

y 1 km J V5 au/ \10 3 Lq/ 
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In other words, the Yarkovsky effect alone can fling an initially 
near-circular asteroid out of the system after about 10 Myr! 
This case however assumes that changes in the asteroid’s rota¬ 
tion state from YORP and close encounters with the star are 
neglected. For highly eccentric orbits, both effects are likely to 
cause a substantial change in the spin state of the body, which - 
in turn - could eventually eliminate the Yarkovsky effect. Also, 
the interplay between the asteroid’s rotation and the stellar wind 
should be investigated. 

In order to compare to the non-Yarkovsky term (which in¬ 
cludes Poynting-Robertson drag and “radiation pressure”), note 
the following standard result (also seen in eauation lll2l) but ap¬ 
plied to GB stars 


( 1 5 AL \ 

1.8 x 10" 5 , 

( P ^ 

\c 2 8nma 2 ) 

Myr ' 

^2 g/cm 3 J 


f R ) 

V a )~ 2 1 

( L \ 

V1 km/ 

V5 aui ' 

\ 10 3 Lq / 
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which is over three orders of magnitude less powerful. 

Although the evolution due to time-dependent radiation is 
more complex than equations lfT08)) -([T09)l suggest (see Appendix 
A), the sheer magnitude of equation (11081) reveals something 
fundamental about the physics in these systems. Possibilities 
include (1) that asteroids in fact escape GB systems en masse, 
with or without the help of mass loss, depending on the stellar 
wind evolution, (2) asteroids survive but are widely dispersed in 
eccentricity and inclination, (3) the temporal variations in the 
individual elements of Q cancel out secular changes in the orbits, 
and (4) the value of k is typically exceptionally low and/or the 
value of Q re f is within a tiny fraction of unity so that Yarkovsky 
is quashed. This degeneracy in interpretation cannot be broken 
without adopting a detailed internal model of the target and in 
particular how its incoming radiation is redistributed and how 
its spin axis changes with time. Such a model is well beyond the 
scope of this paper. 

Further, and unhelpfully, we can pose arguments against all 
four possibilities. Possibility (1) is unlikely because, as described 
in Section 1, asteroids represent the most likely candidate for the 
WD pollution that is currently observed. Possibility (2) might be 
prevented from occurring due to collisions with planetfl, other 
asteroids or smaller bodies, depending on the architecture of the 
exosystem in question. For possibility (3), the equations in Ap¬ 
pendix A demonstrate that even if the off-diagonal nonzero ele¬ 
ments of Q cancel one another, other combinations of elements 
will not cancel. Although the equations were derived assum¬ 
ing position and velocity-independent values of Q, the current 
equations represent an accurate evolutionary picture at system 
snapshots when the seasonal component is neglected. Finally, 
the extreme values of k and Q re f that are required for possibil¬ 
ity (4) do not conform to what we observe in the Solar system. 

We can place these results in the context of the asteroids 
observed in the Solar System. These asteroids, exposed to the 
relatively weak MS luminosity of a IMq star, can experience a 
significant drift in eccentricity (~ 0.1) over 1 Gyr (e.g. Figs. 7-8 


7 Alternatively, distant planet-asteroid interactions due to post-MS 
planet-planet s cattering could repress or excite both eccentricity 
and inclination jVeras &; Armitagell 20051 1 20061 ; I Raymond et al.ll2010l ; 

iMatsumura et al.ll2013l)~ 
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of lVokrouhlickv et al.l2006l 'l . Because these Eros family asteroids 
are located at a ~ 3 au and with R > 3.5 km, this comparison 
yields good agreement with equation (11081) . Note that the semi¬ 
major axis drift in these asteroids over this time period is just 
0.1 au. By estimating the semimajor axis drift from equation 
(fl03l) . we obtain 


1 AL \ ~ 0.81 au / M* \~ 1/2 / p 
c Annina 2 ) Myr \1Mq J \2 g/cm 3 


[ R \ 

rv « r 1/2 i 

( L \ 

y 1 km ) 

\5 au/ 

yio *l @ ) 


( 110 ) 


For a Solar-type star ( L = Lq), this level of drift is within 
the same order of magntiude as the maxi mum that is observed 
and predicted in the Solar system (e.g. iFarinella et all 1 19981 : 
iBottke et al .1 [200(3 : IVokrouhlickv et al]|200fih . 


4-3.5 Averaged equations with no Yarkovsky 


Recall that when Yarkovsky is not active, we are still left 
with (much smaller) radiative perturbations from Poynting- 
Robertson drag and radiation pressure. We compute exactly 
these less-powerful (1/c 2 ) terms to be 


( 
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AL (Qabs + Qref) (2 + 3e 2 ) 
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= 0. 


(114) 


Equations Gnu-dim with a coefficient of 1/c 2 , reproduce 
now-standard results from Wvatt fe Whipple! dl95d) . Their rel¬ 
ative importance is showcased by comparing them to equations 
dH-dH, which instead harbour terms with a coefficient of 1/c. 
Consequently, the Yarkovsky effect when active, dominates the 
target’s secular evolution over other radiative effects. Yarkovsky 
causes changes that are absent from immediately reflected ra¬ 
diation. Equations (I113l) - (lll41) are also important because they 
reveal that when Yarkovsky turns on, a target will drift out of 
the orbital plane that it previously occupied. 


of motions for these forces both during a single orbit and over 
secular timescales (Section 4 + Appendix). 

Several ancillary results are widely applicable to planetary 
systems at all stages of evolution, including pre-MS and be¬ 
ginning MS stages. We summarize these results as (1) a self- 
consistent framework to treat Poynting-Robertson drag, radia¬ 
tion pressure and the Yarkovsky effect ('equation 1291) . (2) an or¬ 
bital element characterization of the Reynolds number for stellar 
winds (equations 1391I49I53I) . (3) the complete set of orbital el¬ 
ement evolution equations for a particle dragged by gas in the 
Epstein and Stokes regimes (equations 1701841) . (4) the complete 
unaveraged set of equations due to radiative perturbations as¬ 
suming position and velocity-independent diurnal and seasonal 
Yarkovsky components (equations 1931971) , and (5) the leading 
order of the averaged set of these equations (IA2IIA6I) when the 
seasonal Yarkovsky component is negligible. This last set, along 
with equations GSED-Gm demonstrate the potential for the 
Yarkovsky effect alone to dominate the eccentricity and inclina¬ 
tion evolution of bodies larger than 1-10 m. 

Detailed models of the time-varying stellar wind and the 
physical properties of an orbiting body are required to determine 
its ultimate fate. Here, we have provided a set of machinery in 
which such models may be incorporated. 
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5 SUMMARY 

Post-main-sequence stellar radiation and winds play a decisive 
role in determining the final orbital states of asteroids, pebbles 
and smaller particles. Any bodies smaller than about 1000 km 
will be significantly affected by forces other than gravity during 
GB evolution. This paper has quantified this main conclusion 
through two avenues: a rough comparison of the instantaneous 
accelerations caused by mass loss, wind drag and radiation (Sec¬ 
tion 3) and a precise detailing of the orbital element equations 
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APPENDIX A: EXPLICIT EXPRESSIONS FOR AVERAGED YARKOVSKY EQUATIONS 

Here we write out the expressions from equations GmMnm, and do so in a form which elucidates how the leading order term (1/c) 
vanishes when the Yarkovsky effect is turned off. We remind the reader that an implicit assumption in these equations is that the 
components of Q are independent of r and v. 

When the Yarkovsky effect plays no role in the motion, then all of the off-diagonal terms of Q vanish, such that Q12 = Q21 = 
Q13 = Q31 = Q23 = Q32 = 0. Further, the diagonal terms become equal, so that Qn = Q22 = Q33. Also, in the subsequent expressions, 
we compress commonly-occurring eccentricity-based quantities as 

k = —2 + e 2 + 2\J 1 — e 2 . (Al) 


The final expressions are 
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where the auxiliary T variables are long explicit functions of (e, We do not write out these functions here but are happy to 

provide them to interested readers. 
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